Chaos in Coulomb-interacting circular billiards 
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We apply a molecular dynamics scheme to analyze chaotic properties of a two-dimensional circular 
billiard system containing two Coulomb-interacting electrons. The interaction strength is varied 
from the noninteracting limit with zero potential energy up to the strong-interacting regime where 
the relative kinetic energy approaches zero. At weak interactions the bouncing maps show jumps 
between quasi-regular orbits, albeit the dynamics is generally chaotic. In the strong-interaction limit 
we find an analytic expression for the bouncing map. Its validity in the general case is assessed by 
comparison with our numerical data. To obtain a quantitative view on the degree of chaoticity, we 
compute the escape rates that, apart from very weak interactions as well as the strong-interaction 
limit, suggest generally chaotic behavior that is independent of the interaction strength. 

PACS numbers: 05.45.Pq,82.40.Bj, 73.21. La 



I. INTRODUCTION 

Classical billiard systems have attracted continuous 
interest for several decades due to their applicability 
to demonstrate chaotic dynamics through (semi) analytic 
and numerical calculations [l|-[3|. On the other hand, 
laboratory experiments on, e.g., microwave billiards 
quantum dots [3], and more recently even graphene 
have rapidly extended the interest in chaos across dif- 
ferent fields in physics. Along this development, billiard 
systems are becoming a key element in the studies of 
classical and quantum chaos both theoretically and ex- 
perimentally. 

Most billiard studies have focused on single-particle 
properties of systems ranging from regular (integrable) 
to chaotic (nonintegrable) systems, including also pseu- 
dointegrable billiards |5j such as regular billiards with 
singular scatterers inside the system. Two-particle bil- 
liards have been studied with hard-sphere contact inter- 
actions in, e.g., rectangular [6| and mushroom-shaped 
cavities. To the best of our knowledge, such studies with 
Coulomb interactions - and with the focus on classical 
chaotic properties - have been restricted to different two- 
dimensional (2D) oscillators (sl-llH. apart from a recent 
molecular-dynamics (MD) study on rectangular billiards 
in magnetic fields [12[. The MD scheme is a computa- 
tionally efficient approach to many-particle billiards that, 
in principle, can be extended to large systems without 
compromising on the numerical complexity of the long- 
range Coulomb interaction. It is noteworthy that the 
Coulomb interaction is a physically meaningful choice 
when considering similar systems in, e.g., quantum-dot 
physics SEIil. 

Here we adopt the MD approach to analyze the classi- 
cal chaoticity of a 2D circular hard- wall billiards with two 
Coulomb-interacting electrons. The system is chosen due 



to its simplicity as a protot3me model for regular billiards 
in the noninteracting Hmit [i| , and the possibility to find 
an approximate expression for the bouncing map in the 
strong-interaction limit. In the intermediate regime the 
system is expected to exhibit chaotic behavior. Due to 
these features the system provides a well-grounded path 
into examinations of both classical and quantum chaos in 
Coulomb-interacting billiard systems. We point out that 
soft billiards are better known in this respect; for exam- 
ple, the two-electron circular harmonic oscillator is reg- 
ular and becomes mixed (partly regular, partl y ch aotic) 
if ellipticity is added in the external potential [11| . 

Our main tool to analyze the chaoticity in an inter- 
acting system is the escape rate. In the noninteracting 
case, Kolmogorov-Arnold-Moser (KAM) islands [1] and 
nonisolated periodic orbits with at least one Lyapunov 
exponent close to zero give rise to a power-law tail in 
the escape rate distribution. On the other hand, fully 
chaotic systems exhibit exponential escape rate distribu- 
tions [ilHigj. 

The paper is organized as follows. In Sec. [Ill we briefly 
introduce the system and our time-propagation scheme. 
In Sec. nil Al we show bouncing maps that demonstrate 
clear signals of chaotic behavior through a large range of 
the interaction strength. At weak interactions, bounc- 
ing maps are found to jump between quasi-regular tra- 
jectories. In Sec. nil Bl we analyze in detail the strong- 
interacting limit and find an approximate expression for 
the bouncing map. The expression agrees with the nu- 
merical result, and at weaker interactions it becomes only 
approximate. Finally, in Sec. IIII CI we assess the degree 
of chaoticity by considering escape rates out of the sys- 
tem. Apart from very weak interactions, we find expo- 
nential escape in a wide range of the interaction strength. 
This indicates strongly uncorrelated trajectories and thus 
chaotic behavior. The paper is summarized in Sec. IV. 
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II. SYSTEM AND METHODOLOGY 

We consider two Coulomb-interacting electrons in a 
circular hard-wall potential. The collisions with the 
boundary are elastic and the system is described by the 
Hamiltonian 



H 



1 



ri -r2 



(1) 



in Hartree atomic units (a.u.) {h = e = rrie = (47reo)~"^ = 
1). Here is the position of the ith electron from the cen- 
ter of the system, and a is a parameter that determines 
the interaction strength. Introducing a in the numerator 
of the interaction term is equal to the scaling the length 
and time by a factor of 1/a and keeping the numerical 
value of the total energy and the velocity the same. In all 
our simulations the total energy of the system is fixed to 
E = 1 and the radius of the circle to = 1/2. The inter- 
action strength is restricted to < a < 1, where a = 
corresponds to noninteracting electrons, and a = 1 corre- 
sponds to electrons being localized at the opposite sides 
of the circle with zero motion. 

To propagate the electrons we use molecular dynamics 
with the velocity Verlet algorithm. The positions 
and velocities of each electron are calculated from 



r{t + At) 
v(t + At/2) 

8i{t + At) 

v{t + At) 



v{t) + U{t)At', 
^F, [r(t + At)]; 



(2) 

(3) 
(4) 



^r{t + At/2) + -a{t + At) At. (5) 



We define 6 and s as the generalized coordinates describ- 
ing the collisions with the boundary. is the angle be- 
tween the velocity vector of the incoming electron and the 
tangent of the boundary, so that < 7r/2 and > 7r/2 
correspond to counterclockwise and clockwise traveling 
directions, respectively. Here s g] — 7r/2, 7r/2] is the ori- 
ented arc length from the chosen origin. 



III. RESULTS 
A. Bouncing maps 

In Fig. [T] we show examples of bouncing maps and elec- 
tron trajectories for a two-electron circular billiard with 
different interaction strengths a. In this case the bounc- 
ing maps consist of 14 000 {a = 10"^) and 5500 {a = 0.2 
and 0.7) collisions with the boundary. The noninteract- 
ing circular system with a = is a well-known exam- 
ple of regular billiards [H represented by straight lines 
in the map (constant bouncing angle) and straight tra- 
jectories forming a star-shaped pattern, where the inner 
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FIG. 1: Examples of bouncing maps (of one of the electrons) 
and trajectories for a two-electron circular billiards with dif- 
ferent interaction strengths: a = 10~^ (up), a = 0.2 (middle), 
a = 0.7 (bottom). Only a small section of the trajectories 
corresponding to the bouncing maps are shown. 



find emerging deviations from this limit as visualized in 
the inset of the upper panel of Fig. [H When the elec- 
trons pass each other the interaction is pronounced and 
we may find "jumps" from one quasi-regular trajectory 
to another one (close- lying parallel lines in the inset). 

In the intermediate-interaction range (middle panel of 
Fig. [T]) the chaoticity of the system is clear, so that the 
bouncing map rapidly becomes completely filled. As ex- 
pected, the distribution of the bouncing map is centered 
at = 7r/2, so that, on the average^ the electrons hit the 
boundary along the normal vector. 

If a is increased above a ^ 0.5 we find that for some 
trajectories the maximum of the probability distribution 
for splits into two. This is visible in the bottom panel 
in Fig. [1] for a = 0.7. However, the splitting is smoothed 
out when a large ensemble of trajectories is taken into ac- 
count. When a is increased further, the system gradually 
becomes (quasi) regular and eventually the bouncing map 
reduces into a one-dimensional curve. In the following we 
carry out analytic calculations in the strong-interaction 
limit. 
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FIG. 2: Two electrons oscillating at the opposite sides of the 
circle when the interaction is strong (a > 0.9). The initial 
positions and velocities are shown. 



B. Strong-interaction limit 
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FIG. 3: Analytic result for the large-a limit [Eq. (|6|)] of the 
bouncing map (solid lines) for a = 0.99, 0.95, and 0.9. The 
corresponding simulated, i.e., numerically exact values are 
shown by points. 



In the srong-interaction limit a ^ 1, the two-electron 
dynamics shows regular characteristics. The electrons 
are confined at opposite sides of the circle as visualized 
in Fig. [21 Here we focus on the special case with total 
angular momentum L = which is conserved due to 
the rotational symmetry. Hence, according to the choice 
of coordinate axes in Fig. [2] we may approximate x = 
ri,x ^ r2,x' We point out, however, that the y coordinate 
does not have mirror-symmetry in the general case. Note 
also the position of s = at x = in Fig. [21 so that 

S ^ [ "^max? "^max] • 

As the second approximation, the electron velocity per- 
pendicular to the edge, can be taken as a constant, 
i.e., it is approximately the same for all possible values of 
s. After straightforward geometrical analysis, taking into 
account the conservation of E and L, we can calculate 
the tangential velocity v\\{s) and further the bouncing 
angle from tan6> = v±/v\\. Thus, we obtain the following 
strong-interaction approximation for the bouncing map 
of electron 1: 



Oi{s) = arctan . 



l2Ko^2Uo-2U,-LlJR^ 
U^-U,{s)^LlJR^ ' 



(6) 



where Kq = Vi^^{0)-\-Vi^y{0) is the initial kinetic energy of 
the electron, Li^o = x{0)vi^y{0) — yi{0)vi^x{^) is its initial 
angular momentum (note that x = xi ^ X2 according to 
our approximation above), Uo = a/[|^i(0)| + 1^2(0) |] is 
the initial potential energy, and R = 1/2 is the radius 
of the circle. Furthermore, Eq. ([6]) has three potential 
energy components corresponding to situations, where 
electron 1 touches the boundary at t = 0, both electrons 
touch the boundary at t = 0, and both electrons touch 
the boundary at s time t. These components, respec- 



tively, have expressions 



a 



|t/2(0)| + v/^?'-x(0)2' 

a 
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a 
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Figure [3] shows the results from Eq. ([6]) for a = 0.99, 0.95, 
and 0.9 (solid lines). The simulated^ i.e., the numerically 
exact values, are shown by points. As a is decreased we 
find gradual deviation from the simulated data. At a = 
0.9 the deviation is already clearly visible. We may thus 
state that Eq. (|6]) provides a reasonable approximation 
for the bouncing map at a > 0.9. This threshold slightly 
depends on the initial conditions; the examples in Fig. [3] 
are chosen such that the deviations between the analytic 
expression and the numerical data are large. 



C. Escape rates 

A natural sign of chaos in a general TV-electron sys- 
tem is the loss of autocorrelation of the trajectory. In 
the present eight-dimensional phase space this is tedious 
to calculate from a single trajectory, and naturally even 
more difficult over all the trajectories in order to per- 
form statistical averaging. However, we can analyze the 
chaoticity of our system indirectly by calculating the es- 
cape rates through a hole on the billiard boundary as a 
function of a. In the single-electron case, systems with 
KAM islands or nonisolated and nonhyperbolic orbits 
have a power-law tail in the escape-rate distributions 
as a function of time or number of collisions with the 
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boundary [TlHTij. In contrast, in fully chaotic systems 
the distribution is exponential. 

We set 10 holes in the boundary covering together 
1/50 of the boundary length - the same fraction as in 
Ref. The escape rates are considered as a function 
of the total number of collisions rather than the propa- 
gation time, since the characteristic time scale strongly 
depends on the interaction strength a. For each a we 
compute 2.5 .. .6 x 10^ respective trajectories with ran- 
dom initial conditions and store the number of boundary 
collisions before the escape. The time steps are chosen 
in the range At = 10~^ . . . 10~^, so that the convergence 
is ensured in every calculation, while the numerical effi- 
ciency is maximized. 

Figure HJ^a) shows the resulting histograms of the 
escape-rate calculations. First, the noninteracting situ- 
ation {a = 0) has a clear power-law tail with P(n) ex 
+ const, where 7 ^ 3.46. In contrast, when 
0.1 < a < 0.5 an excellent fit to the exponential be- 
havior with P{n) = 49^~^/50^ (straight line) can be 
found. This relation results from the system geome- 
try: each collision has the escape probability of 1/50, 
and thus for the n:th collision to lead to escape we find 
P{n) = (49/50)^-^(1/50). This essentially means that 
the correlation with two successive bounces is completely 
lost, and hence the system can be classified as chaotic. 

In Fig. m^b) we have a closer look to the weak- 
interaction limit with a = 10"^. A good fit to the com- 
puted data is obtained with a power-law curve having 
7 ^ 4.08, i.e., the escape is slightly faster than in the 
noninteracting limit. However, at small time scales the 
behavior is very similar to the a = as the trajecto- 
ries essentially follow the same (quasi-) regular patterns. 
These quasi-stable trajectories also give arise to power- 
law tail in the escape-rate distributions for weak interac- 
tions. However, the interaction reduces the lengths of the 
quasi-regular parts of the electron trajectories and thus 
decreases the survival probability (and escape rates) at 
longer time scales. 

Concluding, our results on the escape probabilities 
show that the transition to exponential escape rates - 
and thus to chaotic behavior - is (i) smooth (not abrupt 
as a function of a), (ii) it occurs first at large times (large 
number of collisions) in the histogram, and (iii) it gen- 
erally appears at relatively small values for a. Our tests 
indicate that at a ~ 10"^ the most part of the calcu- 
lated escape-rate histogram is closer to an exponential 
behavior than to the power-law one. We point out that 
the large-a regime is excluded in this analysis due to nu- 
merical reasons: at a > 0.5 we would need to decrease 
the holes due to small-scale motion close to the bound- 
ary, and thus the time step should be decreased as well. 
Hence, for consistency of the results we have focused here 
only on the range < a < 0.5. 
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FIG. 4: (a) Histograms of the escape rates in a two-electron 
circular billiard. The noninteracting case shows power- law 
behavior in the tail (curved solid line), whereas at 0.1 < a < 
0.5 the escape rate is exponential (straight solid line), (b) At 
very weak interactions a — 10~^ we find a mixture of these 
tendencies due to quasi-regular trajectories in the system. 



IV. SUMMARY 



To summarize, we have made a thorough look into 
chaotic dynamics of circular billiards containing two 
Coulomb-interacting electrons with the full range of in- 
teraction strengths (0 < a < 1). Close to both weak- and 
strong-interaction limits the bouncing maps show traces 
of quasi-regular behavior, although the dynamics gener- 
ally appears as chaotic. In the strong-interaction limit 
we are able to find an analytic expression for the bounc- 
ing map that agrees very well with the calculated data 
at a ^ 1. At smaller a the predictive power of the ex- 
pression reduces, although the agreement is reasonable 
down to a ~ 0.9. To assess the chaoticity we have cal- 
culated escape rates as a function of a and found similar 
exponential behavior through a wide range of interaction 
strengths, which indicates chaotic behavior. As the cir- 
cular billiard is completely regular in the noninteracting 
hmit, our finding suggests universally chaotic behavior 
in Coulomb-interacting billiards with hard walls. This 
expected generality is a natural subject of future studies. 
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